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. . , Abstract 

^. 

0\ ■ 

Q^ ■ We consider axisymmetric viscous accretion flows where a fraction / of the 

'""' ! viscously dissipated energy is stored in the accreting gas as entropy and a fraction 

>• I 1 — / is radiated. Assuming a-viscosity we obtain a two-parameter family of 

Q ■ self-similar solutions. Very few such exact self-consistent solutions are known for 

^ ! viscous differentially rotating flows. When the parameter / is small, that is when 

lO ' there is very little advection, our solutions resemble standard thin accretion disks 

'""' . in many respects except that they have a hot tenuous corona above the disk. In 

,__! I the opposite advection- dominated limit, when / ^ 1, the solutions approach 

>■ ■ nearly spherical accretion. The gas is almost at virial temperature, rotates at 

Q^ . much below the Keplerian rate, and the flow is much more akin to Bondi accretion 

than to disk accretion. None of the solutions have funnels. 

We compare our exact self-similar solutions with approximate solutions which 

had been previously obtained using a height-integrated system of equations. We 

find that various dynamical variables such as the radial velocity, angular veloc- 

ity and sound speed estimated from the approximate solutions agree very well 

Q^! with the corresponding spherically averaged quantities in the exact solutions. 

' I We conclude that the height-integration approximation is an excellent one for a 

^ • wide range of accretion conditions, including nearly spherical flows, provided the 

to '. equations are interpreted as spherical averages. 

f^ ' We flnd that the Bernoulli parameter is positive in all our solutions, especially 

>• ■ close to the rotation axis. This effect is produced by viscous transport of energy 

k>< I from small to large radii and from the equator to the poles. In addition, all the 

5^ ■ solutions are convectively unstable and the convection is especially important 

near the rotation axis. For both reasons, we suggest that a bipolar outflow will 

develop along the axis of these flows, fed by material from the surface layers of 

the equatorial inflow. 
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1. Introduction 

In a previous paper (Narayan & Yi 1994, tiereafter NY) we discussed ttie 
potential importance of advection effects in accretion flows. We showed that, both 
at very low and very high optical depths, the energy released through viscous 
stresses in an accretion disk may be trapped within the accreting gas. Most of 
the energy then is advected with the flow as stored entropy. Such advection- 
dominated flows have been found in models of boundary layers in cataclysmic 
variables at low accretion rates (Narayan & Popham 1993) and in models of 
pre-main sequence stars such as the FU Orionis systems at very high accretion 
rates (Popham et al 1993). Advection-dominated conditions may also occur in 
the inner parts of disks around neutron stars and black holes. 

In the analysis presented in NY, we integrated the flow equations in the 
"vertical" direction, i.e. parallel to the rotation axis. Making the usual assump- 
tions of steady state, axisymmetry, and a-viscosity, we obtained a set of ordinary 
differential equations for the gas variables as a function of the cylindrical radius 
R. We showed that these equations have an exact self-similar solution where all 
variables have power-law dependences on R and where the Mach number is inde- 
pendent of R. This solution, which had previously been discovered by Spruit et 
al. (1987), has several interesting and unexpected properties. However, before we 
can explore the consequences of these properties we need first to confirm that the 
self-similar solution itself is real and not just an artifact of the vertical integration 
of the equations. 

Vertical integration is a standard approximation which has been used in 
accretion disk studies from the earliest days. The physical motivation behind this 
approximation is that the vertical thickness of an accretion disk is usually much 
smaller than the local radius, so that the flow velocities are likely to be more or 
less independent of height. It is then reasonable to expect that very little is lost by 
integrating out the vertical coordinate. Unfortunately, the self-similar solutions 
discovered by Spruit et al. (1987) and NY are not thin. The temperature of the 
accreting gas is nearly always close to virial, and the formal vertical thickness is 
comparable to the radius. This inconsistency raises serious questions concerning 
the validity of the solutions and the reliability of the conclusions. 

In this paper, we avoid the height-integration approximation and instead set 
up exact flow equations for steady axisymmetric flow in the r6 plane. The equa- 
tions we obtain are similar to those written down by Begelman & Meier (1982). 
We present here numerical self-similar solutions of the equations. These solu- 
tions are direct generalizations of the height-integrated solutions in NY. Rather 
gratifyingly, we find that all the features of the height-integrated solution are 
reproduced well in the present solutions. In fact, we even find excellent quantita- 
tive agreement between the two approaches, suggesting that the height-integration 
approximation may be much better than previously thought. 

In §2 we write down the equations that we solve along with the boundary 
conditions. Then in §3 we describe numerical results of the equations, discuss 
some of the more interesting properties of the solutions, and make various com- 
parisons with the previously-obtained height-integrated solutions. Finally, in §4 
we discuss the implications of the results, especially for the formation of outflows. 



Appendix A compares our solutions with those obtained by Begelman & Meier 
(1982), and Appendix B discusses the relationship to Bondi (1952) accretion. 



2. Equations of a Steady Axisymetric Advection-Dominated Flow 

We define the isothermal sound speed Cg by 



P= PCs: 



(2.1) 



where p is the pressure and p is the density, and we assume that the gas has a 
fixed ratio of specific heats, 7 = Cp/cy. For convenience we define 
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(2.2) 



We work in spherical polar coordinates r^(/), and write the three components of 
velocity as f^, vg and v^p = Or sin 6*, where O is the angular velocity. The gas 
accretes onto a central mass M and we define the Keplerian angular velocity 

^xir) by 

..M^(^)'". ,2.3) 



Assuming steady state {d/dt = 0) and axisymmetry (d/dcj) = 0), the continuity 
equation gives 
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Integrating this over angle we obtain the net mass accretion rate. 
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We employ the usual ct-prescription for the viscosity, which we write in the 
following form for the kinematic coefficient of viscosity. 
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(2.6) 



where a is a constant. This form is equivalent to the assumption that z/ ~ aCgH 
where H ~ Cs/^k is the "vertical" scale height. Note that z^ is a function of 
position, both because Qk depends on r and because Cs varies from point to 
point. The three components of the momentum equation give (e.g., Mihalas & 
Mihalas 1984) 
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while the energy equation gives 
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(2.10) 

Anticipating the self-similarity form assumed below (eqs 2.11-2.15), we have set 
vq = Q\n the above equations. The left-hand side of equation (2.10) is the gradient 
of the entropy. The right-hand side is the rate of generation of energy through 
viscous dissipation, except that it is multiplied by a parameter /. This parameter 
describes the fraction of the disipated energy which is advected as stored entropy, 
and is therefore a measure of the degree to which the flow is advection-dominated. 
A fraction (1 — /) of the energy is removed through radiative losses. In principle, 
/ could be a function of 6*, but all the results we present here correspond to the 
simplest assumption, viz. / = constant. 

We restrict ourselves in this paper to self-similar flows. We therefore seek a 
solution of the form 



I 



P 

Vr 
V0 



-3/2 



P(^), 



GM 



v{9) =rnK{r)v{9), 



0, 



v^ = rnK{r)n{9), 
Cs = rQ.K{r)cs{9). 
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The form of the solution is obvious. The only lengthscale in the problem is r 
and the only frequency is Qk- Therefore, all velocities must scale with radius as 
rflx- Angular variations at a given radius are modeled through the dimensionless 
functions v{9), n{9) and Cs{0). Given the radial scaling of Vr, the scaling of 
p is uniquely determined by the constancy of M in (2.5), and since r'^pVr is 
independent of r, the continuity equation (2.4) shows that vg = 0. 

By construction, the solution (2.11) - (2.15) automatically satisfies the conti- 
nuity eq. (2.4). Substituting the solution in the momentum and energy equations 
we obtain the following four coupled differential equations in 9: 
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Following NY, we have introduced in eq. (2.19) a quantity e' which we define by 

where / is the parameter we have already introduced in eq. (2.10). 

Equations (2.16) - (2.20) constitute a sixth-order system of ordinary differen- 
tial equations for the four functions v{9), Vt{9), Cs{9), and p{9). The integral (2.5) 
sets the normalization of p{9) and provides one boundary condition. The remain- 
ing boundary conditions are distributed between the equatorial plane, 9 = 7r/2, 
and the rotation axis, ^ = 0. At ^ = 7r/2, we have by symmetry the conditions 

TT dv^dn^dc^^dp^ 

2 d9 d9 d9 d9 ^ ^ 

At ^ = 0, we insist that the solutions be well-behaved and non-singular. This 
leads to the conditions 

The last condition follows from eq. (2.19). Not aU of the conditions (2.21), (2.22) 
are independent. We choose a convenient subset of these conditions and solve the 
differential equations (2.16)-(2.19) using a numerical relaxation technique (e.g.. 
Press et al. 1992). 

Under certain conditions, the sixth-order set of equations (2.16)-(2.19) re- 
duces to a second-order system. We discuss this simplification in Appendix A 



and compare our work to a previous analysis of this problem by Begelman & 
Meier (1982) 

Technically, equation (2.19) allows two possible boundary conditions on v at 
^ = 0, viz. V = which is the one we have given in eq (2.22) and v = —e' /2a. 
We have tried to obtain rotating solutions which satisfy the second condition 
and have been unable to find any. However, the boundary condition v = — e'/2Q! 
does allow a purely spherical non-rotating inflow solution which we discuss in 
Appendix B. This solution is a generalization of Bondi flow to the case when 
there is viscosity. 



3. Results 

3.1 Typical Solutions 

We have obtained numerical solutions of eqs. (2.16) - (2.19) for a variety 
of values of the viscosity parameter a and the thermodynamic parameter e' (de- 
fined in eqs. (2.6) and (2.20)). Figure 1 shows a typical sequence of solutions 
corresponding to a = 0.1 and e' = 0.1, 1, 10. These solutions may be con- 
sidered either as a sequence of fully advection-dominated flows (/ = 1) with 
7 = 1.6061, 1.3333, 1.0606, or as flows with a fixed value of 7 and with a 
sequence of decreasing / or increasing cooling. 

The four panels in Fig. 1 show the variation with polar angle 9 of various dy- 
namical quantities in the solutions. The top left panel displays the dimensionless 
angular velocity Q{9). Rather surprisingly, we find that O is nearly independent 
of 9 in each solution, varying by only ~ 10% from ^ = to ^ = 7r/2. Radial shells 
therefore rotate more-or-less rigidly, but of course there is differential rotation 
between neighboring shells. The actual value of O varies significantly from one 
solution to another, changing from O ~ 0.2 at e' = 0.1 to Q ~ 0.9 at e' = 10. 
The scahng of O with e' follows eq (3.2). Note that Q oc (e')^/^ for e' < 1. This 
implies that O ^ as 7 — > 5/3. 

The top right panel shows the radial velocity profiles v{9) of the solutions. 
The velocity is zero at ^ = (this is a boundary condition) and maximum at 
9 — n/2. We find that v is essentially independent of e' for e' <^ 1 and varies as 
V oc 1/Ve' for e' ^ 1 (see eq 3.1). 

The bottom left panel shows profiles of the density p{9). In the e' = 0.1 
solution p varies by only ~ 10% from ^ = to ^ = 7r/2. The solution therefore 
corresponds to a nearly spherical configuration. This is demonstrated in Figure 
2 (top left) where we display isodensity contours in the meridional plane. The 
resemblance of this solution to a star is striking, but of course it is not a normal 
star, since it involves a steady accretion flow. In any case, it is quite clear that this 
solution is very definitely not a "disk" in the usual sense. The density contrast 
between p(0) and p(7r/2) increases with increasing e', becoming a factor ~ 2 at 
e' = 1, and a factor ~ 50 at e' = 10. The isodensity contours of these solutions 
are shown in the top right and bottom left panels of Fig. 2. The e' = 1 solution 
looks like a rotationally flattened star, while the e' = 10 solution is beginning to 
resemble a standard thin disk. A value of e' = 10 normally implies a small value 
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of / (unless 7 is close to 1), and so the solution represents a case where there is 
significant cooling. This is precisely the limit where we expect the fiow to occur 
in a thin disk. 

The v{9) and p{9) profiles both peak at 6* = 7r/2. Therefore, in all our 
solutions the bulk of the accretion occurs along the equatorial plane, and the 
accretion rate goes to zero along the rotation pole. 

The bottom right panel of Fig. 1 shows the variation of c^, or equivalently 
the gas temperature, with 9. In the e' = 0.1 solution, c^ is almost independent 
of 6*, and the pressure p = pc^ too is independent of 9. In this solution, the ro- 
tation is highly sub-Keplerian and so hydrostatic equilibrium requires primarily 
a balance between gravity and the pressue gradient. Since gravity acts in the 
radial direction, the pressure gradient too is almost radial. The e' = 1 and espe- 
cially the e' = 10 solutions have larger temperature and pressure variations with 
9 and have non-radial pressure gradients. This is to be expected given the more 
rapid rotation of these solutions and the increasing importance of centrifugal 
acceleration. 

An interesting feature of the large e' solutions is worth emphasizing. As 
already mentioned, these solutions have efficient cooling and therefore resemble 
thin disks. Nevertheless, in all cases there is a low density corona above the disk 
which is at nearly virial temperature. This is illustrated by the e' = 10 solution 
in Fig. 1. The hottest temperature is achieved at the rotation poles, 6* = 0, n. 

In addition to the above examples with a = 0.1, we have calculated a number 
of solutions with other values of a. For a <til, v scales as a, but except for this, 
solutions with the same e' but with different values of a are virtually indistin- 
guishable from one another. There are more significant variations when a exceeds 
unity. However, such large values of a are probably unlikely (e.g. Narayan, Loeb 
& Kumar 1994, Hawley, Gammie & Balbus 1994), and we have not explored this 
region of parameter space. 

In the advection-dominated limit, where e' ^ 1, our solutions have very 
simple scalings, viz. v ~ —a, O ~ (e')^'^, Cg ~ 1. These scalings arise as follows. 
Since the cooling is inefficient, all the viscous energy is stored in the accreting 
gas, and this means that the thermal velocities approach virial speeds, i.e Cs ~ 1. 
Comparing the various terms in eq (2.18) we see that v has to scale as —ac^. This 
gives V ~ —a for a virial gas. The physical reason for the scaling is that the radial 
velocity in an accretion fiow is determined primarily by the rate at which angular 
momentum is removed from the gas and this depends on the viscosity coefficient u. 
In the a prescription (eq 2.6), we have z/ ~ a in scaled units and this therefore 
implies v ~ —a. Finally, the scaling of O arises through the energy equation 
(2.19). The left-hand side of this equation is the product of the radial gradient 
of the entropy and the radial velocity, and therefore represents the steady state 
rate of change of entropy of a parcel of accreting gas. Recall that our solutions are 
by construction self-similar with p oc r"^'^ and p oc r~^'^. If 7 is exactly equal 
to 5/3, then a fiow with these radial dependences is automatically isentropic. 
However, for 7 7^ 5/3, there is an entropy gradient in the fiow such that the 
entropy increases inwards whenever 7 < 5/3, i.e. e' > 0. This means that the 
entropy of each accreting gas element increases with time at a rate proportional 
to e'. The entropy has to be generated of course by the viscous energy dissipation, 
described by the four terms in the right of eq (2.19). Usually, the dissipation is 



dominated by the two terms proportional to fl'^ amd (dfl/dO)'^ which arise from 
the r(j) and (pz components of the shear stress. Thus the dissipation is proportional 
to O^. In order to achieve self-consistency, the magnitude of O in the self-similar 
solution is adjusted such that the dissipated energy exactly matches the energy 
that is required to maintain the self-similar entropy gradient. Since the latter is 
proportional to e' and since v ~ —a, we thus see from (2.19) that we require 
O ~ (e')^'^. This is exactly what we see in our numerical solutions. 

A somewhat interesting feature of eq (2.19) is that even when O ^ 0, 
the right hand side still remains non-zero. This is because of the terms v"^ and 
(dv/dO)'^, the first of which is the viscous dissipation arising purely from the ge- 
ometric convergence of the flow due to the spherical geometry, while the second 
represents the rz component of the shear stress. These two terms imply a certain 
minimum level of viscous dissipation even in a very slowly rotating flow. This 
dissipation will of course cause the entropy to increase inwards, and by the argu- 
ments given above we see that e' has to be greater than a certain minimum value. 
Indeed, we have discovered from our numerical experiments that self-similar so- 
lutions exist only for e' > Ca^, where C is a constant of order unity. The limit 
is exactly of the form we expect from eq (2.19). Thus for large values of a, the 
parameter 7 needs to deviate significantly from 5/3 and/or / needs to be quite 
different from unity in order to have a self-similar flow. On a related point, we 
note that there is a second branch of solutions which corresponds to non-rotating 
purely spherical accretion. These solutions are quite distinct from the O — i> 
limit of the rotating solutions discussed here and are in fact closely related to 
Bondi (1952) spherical accretion. We discuss these solutions in Appendix B. 

In addition to the solutions described so far, we have found (for somewhat 
large values of e'/a^) other solutions where O reverses sign one or more times 
as a function of 9. These higher-order solutions come in two parities. Solutions 
with even parity have an even number of nodes in Q{9) between ^ = and 
9 = TT. These solutions satisfy the boundary conditions (2.22) at the equator, 
and have rotation proflles which are symmetric between the two hemispheres. 
The lower right panel of Fig. 2 shows isodensity contours of one such solution 
with two nodes where the flow down the two poles rotates in one sense while 
the flow in the equatorial plane rotates with the opposite sense. Solutions with 
odd parity have an odd number of nodes, and the rotation proflles in the two 
hemispheres are reversed with respect to each other. These solutions satisfy the 
boundary condition O = at 6* = 7r/2. We do not expect any of these higher- 
order solutions to be relevant except in rare cases where the initially infalling gas 
happens to have reversals in the sign of the angular momentum as a function 
of 9. In the rest of the paper we restrict ourselves to the nodeless fundamental 
solutions. 

The results described so far correspond to the particular form of viscosity 
given in eq. (2.6). To check how sensitive the results are to the viscosity pre- 
scription, we obtained solutions corresponding to a second law. By dimensional 
analysis we see that self-similarity is possible only if u scales with radius as r^'^. 
Therefore, as our alternate prescription we used u = aCsV. We found that the 
solutions with this viscosity law are very similar to those described above. There- 
fore, none of our results are special to the particular viscosity prescription we 
have adopted. 
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3.2. Comparison with Height-Integrated Solutions 



One of the primary aims of this study was to check the validity of the 
height-integrated approximation for advection-dominated flows. As we have seen, 
accretion flows become very nearly spherical when they are advection-dominated 
(/ -^ 1). From this it might appear that the height-integrated approximation, 
which is based on a disk-like picture of the flow, would be particularly inappro- 
priate in this limit. We investigate the issue quantitatively. 

Spruit et al. (1987) and NY showed that the height-integrated equations 
have the following analytical self-similar solution. 
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where the second relation in each equation refers to the limit a <^ 1, and the 
subscript h is to remind us that these expressions correspond to the height- 
integrated approximation. The function g{a, e') is given by 
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We now take the exact solutions calculated in this paper and extract from 
them three fiducial values each of f, O and c^. First, we consider the values of 
the variables at the mid-plane, 9 = 7r/2; we refer to these values as {v)m, (^)m 
and {c^)m- As our second estimate, we compute spherically-averaged values, e.g. 
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with {Q)g and {0^)0 defined similarly. Finally, we compute ^-averaged values, 
which correspond to cylindrical averages parallel to the rotation axis. In this 
case, we define the averages according to 
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(3.6) 



(3.7) 
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Figure 3 compares these three fiducial values of f, O and c^ from the exact 
solutions with the corresponding height-integrated values, (f)^, (fi)^, {c'^)h, for 
a range of e' extending from e' = 10~^ to 10^. The results correspond to small 
values of q: < 0.1; in this regime, O and c^ are independent of a while v is simply 
proportional to a. 

When e' ^ 1, we are in the cooling-dominated regime where the fiow resem- 
bles standard thin disk accretion. In this limit, we expect height-integration to be 
quite accurate, and indeed we do find that the height-integrated estimates {v)h, 
{fl)h, {ci)h agree well with all three estimates obtained from the exact numeri- 
cal solutions. The only discrepancy is between {c'^)h and (Cg)m, but this can be 
understood. These solutions have coronae as discussed in §3.1, and naturally the 
midplane value of c^ is smaller than either the spherical average or the ^-average. 

When e' falls below unity the fiow becomes quite spherical, and we would ex- 
pect height-integration to be less valid. As expected, we find that the 2-averaged 
values of V, Q and c^ differ significantly from the height-integrated estimates. 
For instance, at e' = lO'^, we find {v)^ = 0.32{v)h, {^)z = 0.42(0)/,, (c^)^ = 
0.46(Cg)/i. Thus, height-integration leads to fairly large errors at the level of fac- 
tors ~ 2-3. 

However, Fig. 3 reveals a surprise, viz. that the height-integrated estimates 
agree very well with both the midplane values and the spherically-averaged values 
of the exact solutions. For instance, at e' = 10~^, we find {v)m = l-28(t')/,, 
(O)^ = 1.13(0)/i, {c^)m = l-00(Cg)/i for the midplane values, and {v)0 = 0.82{v)h, 
{'D,)g = 1.09(O)/i, {c^)9 = 1.00(Cg)/i for the spherical averages. The spherical 
averages show particularly good agreement with the height-integrated estimates 
over the entire range of e' from to infinity, with no error exceeding 20%. 

This result suggests that one should interpret the height-integrated equations, 
not as averages over cylindrical height z, but rather as averages over spherical 
polar angle 6* at a fixed r. Once this is done, height-integration is a good approxi- 
mation even for nearly spherical fiows. The height-integrated equations, especially 
in the form of the so-called "slim disk" equations (Abramowicz et al. 1988), have 
become popular in recent years for modeling the dynamics of accretion fiows. Un- 
til now, it has not been clear exactly how slim a disk has to be in order for the 
equations to be valid, and also exactly what kind of an average the solutions 
represent. Based on the results presented here, we suggest that the slim disk 
equations may be applied virtually to any accretion fiow around a point mass, 
however non-slim the fiow may be, and that the results should be interpreted as 
spherical averages rather than as ^-averages. We caution, however, that this sug- 
gestion is based on the properties of a very special class of self-similar solutions, 
and needs to be tested on non- self- similar fiows. 

3.3. The Bernoulli Parameter 

The Bernoulli parameter Be, defined as the sum of the kinetic energy, the 
potential energy and the enthalpy of the accreting gas, is of interest in accretion 
fiows because it measures the likelihood that outfiows or winds may originate 
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spontaneously (NY). An adiabatic flow has a constant Be along streamlines. If 
Be is positive for any of the accreting gas, then this gas can potentially reach 
infinity with a net positive kinetic energy. 

Let us normalize Be in our self-similar solutions by the square of the local 
free-fall speed, Vff = ^k^, and consider the dimensionless parameter 

-- ^^ ^y^ + hnsmef-1 + ^cl (3.9) 



nj^r^ 2 2' ' 7-1 

The parameter 6 is a function of Q but independent of r. In NY we showed 
that the Bernoulli parameter is positive in height-integrated advection-dominated 
flows, and suggested that this may explain the frequent occurrence of outflows 
and winds in many accretion systems. However, because the result was obtained 
through the height-integrated equations, it had to be treated with caution. We 
now consider the behavior of h in the exact self-similar solutions obtained in this 
paper. 

Fig. 4 shows 6(^) for a sequence of solutions. All the solutions have the 
same e = 0.333 (7 = 1.5), but the advection-domination parameter / varies over 
the range / = 1, 0.446, 0.33, 0.033, 0.0033, 0.0011, which means that e' = 
0.33, 0.75, 1, 10, 100, 300. We confirm the basic result of NY that advection- 
dominated flows have positive h. Furthermore, we see that for / > 0.446, 6(^) > 
at all Q. This means that for flows that are so highly advection-dominated, the 
entire gas has positive Bernoulli parameter. Even in the limit of extreme cooling 
(/ ^ 0), there is always some gas near the rotation axis [Q -^ 0) which has a 
positive h. The reason for this is the corona which we mentioned in §3.1. Even 
though the disk cools efficiently and is geometrically thin, it still has a hot virial 
corona above it and this region of the fiow can acquire a positive h especially 
close to the rotation axis. 

These results make the connection between accretion and outflows much 
stronger than we suspected based only on the height-integrated work of NY. 
A case can be made now that perhaps all accretion flows, whether advection- 
dominated or cooling-dominated, are capable of producing outflows. The main 
difference between the two kinds of flow may be only a quantitative one, viz. the 
former probably produce more powerful outflows than the latter. Also, it would 
appear that outflows will prefer to form along the rotation axis, since this is 
where h is most positive in all cases. A bipolar morphology is thus natural. The 
argument is made even stronger by the results on convection discussed in the 
next section. 

It is important to emphasize that the positivity of h does not imply a violation 
of energy conservation. In viscous accretion flows, the energy content of a parcel 
of gas is modifled by viscous transport of energy from one radius to another and 
from one Q to another. Let us consider only the radial flux for simplicity. The 
radial flux of energy transported by viscosity is equal to the quantity (the radial 
angular momentum flux) x (the angular velocity) ~ ^'^%- This energy flux is 
directed outward and has a negative divergence. Therefore, it deposits energy at 
each radius. For a given M, the rate of deposition of energy is maximum in the 
case of a standard thin accretion disk, where v^ = Q.kt. However, because these 
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disks also cool efficiently, the deposited energy does not produce an enhancement 
of 6 (in fact, b = —0.5), but rather leads to an enhancement of the locally radiated 
flux. This explains the well-known result that the radiative flux emerging from 
a thin disk (except near the inner edge) is three times larger than the rate at 
which gravitational energy is released locally (e.g. Frank, King & Raine 1989). As 
the parameter / increases and the flow becomes more advection-dominated, the 
viscous energy flux actually decreases in magnitude because V(f, becomes smaller 
than 0/^r. However, because the cooling is also less efficient, a larger fraction of 
the energy is retained by the gas and therefore, paradoxically, b actually becomes 
larger. As we have seen, when / exceeds a critical value, b can actually become 
positive over the entire flow. 

The 9 component of the viscous stress causes angular redistribution of the 
energy at a given r. The direction of this flux is such as to enhance b at the 
poles relative to the equator. The result is that the pole always has a positive b 
even in the limit of small /. 

The self-similar solution is very special in that it has an inflnite source of 
energy at r = which can be transported outward by viscosity. This explains 
how the entire solution can have a positive b. If we consider a non-inflnite flow 
which is terminated at a flnite inner radius r^, then for radii close to r^ the viscous 
energy flux will be less than the self-similar value. The gas near the inner edge 
will therefore have negative b. However, once we are reasonably far from the inner 
edge (say r > few x n), the self-similar value of b will be achieved and beyond 
this radius the flow will be indistinguishable from the self-similar form. Overall, 
the deflcit of b near the inner edge will compensate for the positive b elsewhere, 
ensuring that energy is conserved globally. 

3.4. Convection 

In NY, we showed that advection-dominated flows have entropy increasing 
inwards and therefore that these flows are intrinsically unstable to convective in- 
stabilities. This point was made earlier by Begelman & Meier (1982). We discuss 
here the role of convection in the solutions described in this paper. 

The left-hand side of eq (2.19) is proportional to vTds/dr, where ds/dr is 
the radial derivative of the entropy. Since v is negative and the right-hand side 
of eq (2.19) is positive (it consists only of dissipation terms), we see that ds/dr is 
negative at all 9, and therefore that the flow is convectively unstable at all angles. 
In any simple theory of convection we expect the convective flux to be parallel to 
the local pressure gradient since this is the direction associated with the buoyancy 
force. (In principle, the flux could be in a different direction in a rotating flow 
because of anisotropic transport, but we ignore this complication.) We have seen 
earlier that whenever the flow is advection-dominated (/ ^ 1, e' < 1) the flow is 
nearly spherical and the pressure gradient is almost radial. Therefore, we expect 
convection to act primarily in the radial direction. This is very different from 
the vertical convection which is usually discussed in the context of thin accretion 
disks (e.g. Ryu & Goodman 1992). 
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When there is convection we can expect it to have a back reaction on the 
basic flow. To estimate the magnitude of the effect, let us rewrite the energy 
equation (2.19) with the convective contribution included: 



3e pvrcj , fapcj 



+ ^ ^ + - sm^ en^ + sm^ { —-] 



r'^ r^ \d9 J 4 \d9 J 

(3.10) 
Here Fc is the convective energy flux, and p, Vr, fi, Cg refer to the physical 
variables and not the scaled functions deflned in eqs (2.11)-(2.15). The term on 
the left of eq (3.10) is the divergence of the advected entropy flux, pTv ■ Vs, 
which we refer to as the "advection term." The flrst term on the right is the net 
deposition of energy due to the inflow of convective energy flux into the gas — 
the "convection term" — and the second term is the energy deposition due to 
viscous dissipation (reduced by our usual factor /). 

To estimate the magnitude of the convection term, we note that we expect 
the convective flux to be proportional to the entropy gradient with some effective 
diffusion constant Kc- Further, we may write Kc approximately in a form similar 
to the ct-prescription for viscosity. Thus, we have 

2 

Fe ^ -K^pT{Vs ■ r)r ^ -^^pT{Vs ■ f)r, (3.11) 

where we have assumed that Fc is parallel to the unit radial vecotr f because 
the pressure gradient is primarily radial. Since Vs ■ f is negative we see that 
the convective flux flows outwards. In an advection-dominated self-similar flow 
Cg < 20|^r^/5 (see eq 3.3), and Fc is proportional to r~^ . Therefore, the con- 
vection term becomes 

-V ■ Fe = ^ < la^epclVLK. (3.12) 

r 5 

We see that this term has the same form as the advection term on the left of 
equation (3.10). We can therefore compare the two terms directly. Figure 5 shows 
the ratio of the convection and advection terms for one of our advection-dominated 
solutions, assuming ac = a/2. The shape of the curve may be understood as 
follows. For small e', the radial velocity v is given quite accurately by 

v{e)^-^asm'^enKr. (3.14) 



Therefore, we expect 



4 
Convection 8 a 



Advection 15 sin 6* a 



(3.15) 



The ratio is smallest in the equatorial plane {9 -^ 7r/2) and diverges towards the 
rotation axis {9 ^ 0, tt). 

On general grounds, we would argue that ac cannot exceed a. The reason is 
that a, which describes the viscous transport of angular momentum, is expected 
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to have contributions from many sources of viscosity such as magnetic stresses, 
convection, and other fluid instabihties. Convective energy transport on the other 
hand depends only on the strength of the convection. We therefore expect < 
otc < ct, with the upper equality being achieved only in the limit of massively 
efficient convection with no source of viscous stress other than convection. In 
fact, if magnetic stresses always saturate at equipartition (e.g. Gammie, Hawley 
& Balbus 1994), then the inequality becomes < etc < a/2. 

Figure 5 and eq (3.15) show that for a wide range of ^, the convection term is 
smaller than the advection term, so that convection has only a moderate effect on 
the flow. To see this in more detail, we take the convection term to the left-hand 
side of eq (3.10) and combine it with the advection term (NY). This gives 

... „ . 3e pvci / 8 ac\ ,^ ^^. 

Advection — Convection = — I 1 ^^ . (3.16) 



15sin^6' a 

With this modification, eq (2.19) continues to be valid except that we have to 
rewrite e' as 

e' = -,U-^^-)- (3.17) 

The new factor in parentheses modifies the value of e'. So long as this factor is 
positive, it has only a minor effect on the solution. The value of e' is reduced (in 
an angle-dependent way because of the sin 9) and the fiow becomes effectively 
more advection-dominated. The overall structure of the fiow, however, remains 
basically unchanged. Since the convection term is smaller than the advection 
term, the timescale on which convection changes the entropy profile is longer 
than the advection timescale, and there is just not enough time for convection 
to have a large effect. 

There is however a critical angle Ocriti whose value is given by eq (3.16), 



Or.rit ~ sin 



c \ 1/2' 



15a 



(3.18) 



below which convection dominates over advection. This means that for 9 < Ocrit 
the fiow is strongly modified. Indeed, at these angles, convection will overwhelm 
the advection and the whole self-similar fiow will break down. At these angles, we 
expect that convection will rapidly transfer entropy outwards. Since the Bernoulli 
constant is positive at these angles (see §3.3) the gas at large radius will acquire 
a great deal of positive energy and in all likelihood will fiow out supersonically. 
We thus have a plausible scenario for the formation of bipolar outfiows. Our 
qualitative picture is that the accretion occurs primarily in a thick equatorial belt 
with Ocrit < < t: — Ocrit- Over this zone the solution is almost of the self-similar 
form described in this paper, and although it is convective the convection has 
only a minor effect on the parameters of the fiow. However, for < Ocrit and 
> TT — Ocrit, the convection becomes strong enough to transfer energy outward 
and to disrupt the self-similar accretion. This gas has a positive Bernoulli constant 
(§3.3), and therefore very likely will be driven out in an outfiow. Of course, we 
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still don't have a self-consistent description of the outflow region, since all we have 
shown is that our self-similar advection-dominated solution is violently unstable 
in the region around the rotation axis. Nevertheless, we feel that our picture is 
plausible and that some outflows at least may well be generated in this manner. 

Figure 5 shows that for ac = a/2, we have 9crit ~ 33°. The outflow region 
is therefore not very wide. Moreover, the amount of mass involved in the outflow 
region is quite small. Let us assume that the mass outflow rate Mout is equal to 
the mass accretion rate in the original self-similar solution for 9 < Ocrit- While 
the density p is essentially independent of 9 in an advection-dominated flow, the 
accretion velocity varies as t; oc sin 9 (eq 3.14) and is very small in the region 
around the axis. From our numerical solution we calculate that Mout ~ 0.017M 
for ac = a/2. Therefore, only a small fraction of the accreting mass participates 
in the outflow. Indeed this estimate of Mout is probably too high since it assumes 
that the original self-similar flow has been set up in the polar regions and is then 
turned round into an outflow. In practice we imagine that the outflow will clear 
out two conical regions and that the mass for the outflow will be supplied by 
pressure gradients in the upper layers of the equatorial inflow. The mass outflow 
rate will then be even smaller than our estimate. 

In any case, the important point is that, in the very nature of the flow, 
the mass in the outflow region acquires the most positive energy (or Bernoulli 
constant) at the expense of the rest of the accreting material. It is therefore 
ejected with a speed comparable to the free fall speed at the radius from which 
the outflow originates. The most energetic ejected material will have a speed at 
inflnity comparable to the virial speed at the surface of the accreting star. 



4. Summary and Discussion 

The main aim of this investigation was to obtain axisymmetric self-similar 
advection-dominated flow solutions in three dimensions without using the height- 
averaging approximation. We have succeeded in this enterprise. The solutions we 
have presented are obtained by numerically solving a general sixth-order system 
of differential equations, where the only serious approximation we have made is 
the use of an isotropic a viscosity. To our knowledge, these solutions are one of 
the very few fully self-consistent, axisymmetric, rotating steady state flows known 
with non-trivial viscous interactions. Related previous work has been published 
by Begelman & Meier (1982), Liang (1988) and Henriksen & Valls-Gabaud (1994). 
However, some of these other solutions have unphysical boundaries where the as- 
sumptions break down, whereas our solutions describe fully consistent equilibrium 
flows which flU the entire r9 plane. 

Our solutions span a two-parameter family labeled by the viscosity parameter 
a (see eq. 2.6 for the deflnition) and a thermodynamic parameter e' (eq. 2.20). 
The latter is a function of the ratio of speciflc heats 7 of the accreting gas and the 
fraction / of the dissipated energy which is advected with the gas. For a given 
a, solutions with large values of e' behave like standard thin disks, as might 
be expected since these solutions correspond to / — > and so advect very little 
energy. In the opposite advection-dominated limit, which corresponds to / ^ 1 
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or e' < 1, our solutions describe nearly spherical flows which rotate at much 
below the Keplerian rate. These advection-dominated solutions have very similar 
properties to the approximate solutions derived in NY. 

We emphasize that the solutions we find are not the usual tori with steep 
funnels that are normally considered when "thick accretion disks" are discussed 
(e.g. Begelman & Meier 1982, Frank et al. 1992). Our solutions are much more 
like slowly rotating stars, except that they are not static but involve viscously- 
driven settling flows. The question now arises: do our equations permit solutions 
with empty funnels, as described by Begelman & Meier (1982)? As we discuss 
in the Appendix, this must be considered an open question at this time. 

We have compared the numerical three-dimensional solutions of this paper 
with the approximate two-dimensional height-integrated solutions described by 
Spruit et al. (1987) and NY. Because some of our solutions are nearly spherical, 
we might expect the height-integrated solutions to be in error by quite large fac- 
tors. Instead, we flnd that the height-integrated solutions agree very well (errors 
< 20%) with spherically-averaged quantities from the exact numerical solutions. 
This is very encouraging since it means that simple height-integrated equations 
such as the "slim disk" equations (Abramowicz et al. 1988) may be valid over a 
much wider range of conditions than suspected before. The result also clarifles 
how to interpret the height-integrated equations. One should consider these equa- 
tions to represent the properties of spherically- averaged quantities in the accretion 
flow (averages over spherical polar angle 9 at fixed spherical radius r) rather than 
of vertically averaged quantities (averages over vertical height z at fixed cylin- 
drical radius R). This interpretation is, however, based on a very special class of 
solutions and should be checked on less special fiows. 

Advection-dominated accretion fiows have a very unique feature in that they 
are characterized by positive values of the Bernoulli parameter b (eq 3.9). This re- 
sult was discovered by NY using the height-integrated equations, and is confirmed 
here with more exact calculations. We interpret a positive Bernoulli parameter 
as an indication that the accreting gas may be able spontaneously to generate 
winds and outfiows. In some of our highly advection-dominated solutions we find 
that the entire accreting gas has a positive b (e.g. the solutions with / > 0.446 in 
Fig. 4). These fiows may be particularly susceptible to violent outfiows. Rather 
surprisingly, even solutions with efiicient cooling (/ -^ 0) have hot low density 
coronae with positive b. This might mean that even regular thin accretion disks 
may always be accompanied by advection-dominated coronae which can drive low- 
density winds. Further, both in the advection-dominated and cooling-dominated 
limits, b is maximally positive along the rotation axis. This suggests that when- 
ever there is an outfiow, it is likely to have a generic bipolar morphology. 

Another interesting property of advection-dominated fiows is that they are al- 
ways convectively unstable, as was noted originally by Begelman & Meier (1982). 
Because advection-dominated fiows have almost spherical isobars, we argue that 
the convection will occur almost purely in the radial direction. This is very dif- 
ferent from the "vertical convection" which has been studied in the context of 
thin accretion disks (e.g. Ryu & Goodman 1992). 

Convection in advection-dominated fiows may be a source of turbulent vis- 
cosity. Although viscosity is required in order to have any accretion at all, for 
long it was unclear what the source of the viscosity in accretion disks may be 

16 



since no linear hydrodynamic instability could be identified in Keplerian disks. 
With the recent work of Balbus & Hawley (1991), interest therefore shifted to 
magnetic stresses generated through an MHD instability. The only model that 
we are aware of that explains accretion disk viscosity through a purely hydrody- 
namic model is that of DubruUe (1992) who uses a finite amplitude instability. 
This model, however, results in a very small a <^ 1. The convective instability 
which we identify in our advection-dominated flows is a true linear instability. 
Furthermore, if the flow is self-similar, then the convection also will automatically 
be self-similar, and this means that the convective viscosity will necessarily have a 
form similar to the a-scaling given in eq. (2.6). Our models are thus rather close 
to being perfectly self-consistent in the sense that the viscosity is not imposed in 
an ad hoc manner but could, if necessary, itself be determined self-consistently 
from an instability in the flow. 

By comparing the convective time scale with the advective time scale, we 
flnd that for most regions of our solutions, convection acts only as a moderate 
perturbation which changes the parameters of the flow by a modest amount 
but does not destroy the overall structure of the flow. However, for a range of 
angles around the rotation axis, 6 < Ocrit ~ 30°, the convection is capable of 
overwhelming the advection and may transport entropy outwards more rapidly 
than it can be carried in by accretion. Since this region of the flow is also the 
zone with the largest Bernoulli parameter, we suggest that the violent convection 
and the positive h will together act to reverse the local accretion flow into a 
bipolar outflow. The speed of the ejected material at inflnity will be of order the 
free-fall speed at the radius of origination of the material. We imagine that the 
outflow will be fed by surface material from the equatorial inflowing gas. 

We feel that this scenario is a plausible and generic mechanism to produce 
outflows. For the mechanism to work, we require the outflowing material to re- 
main adiabatic and to retain its Bernoulli parameter long enough to be accelerated 
beyond the escape velocity. Given this requirement, we suggest that an outflow 
can be created merely with viscous and convective redistribution of energy in 
the manner we have described, without any additional agencies. Hydrodynamic 
models of outflows have been discussed before by Eggum, Coroniti & Katz (1988), 
Liang (1988) and Henriksen & Valls-Gabaud (1994). In the former two papers, 
the outflow is accelerated by radiation pressure, while magnetic stresses appear 
to play a role in the last work. Radiation and magnetic flelds doubtless help, but 
we suggest that outflows are possible even without them. 

Note that we do not at this time have a truly self-consistent description of 
the outflow. By insisting on self-similarity and a constant M at all radii, we are 
forced to set vq = 0, and because of this our solutions do not have the option of 
diverting any of the accreting mass into an outflow. The condition of constant 
M is relaxed by some other authors (e.g. Liang 1988, Henriksen & Valls-Gabaud 
1994), and they do flnd solutions with outflows. It is not clear, however, whether 
their outflows are driven by the same mechanism that we suggest. It would be of 
interest to reformulate our problem with a non-constant M to investigate whether 
or not self-similar outflows are possible with these equations. 

The principal characteristic of an advection-dominated accretion flow is that 
most of the dissipated energy is stored within the gas. This can happen at 
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very low accretion rates when the gas density is so low that cooling through 
brenisstrahlung is inhibited. Such an effect was seen in detailed numerical models 
of boundary layers in cataclysmic variables (Narayan & Popham 1993). Alterna- 
tively, advection-domination may occur when the accretion rate is extremely high. 
In this case, the optical depth of the flow becomes so large that the radiation is 
unable to escape in less than a flow time. Popham et al. (1993) found that this 
may happen in FU Orionis systems. Advection-domination has been considered 
also in disks around supermassive black holes. Rees et al. (1982) discussed ion 
tori at very low accretion rates, while Begelman (1978) discussed radiation trap- 
ping at very high M where the photon diffusion time scale becomes longer than 
the radial inflow time scale (or advection time scale). 

If advection-dominated conditions are present over a sufficiently wide range 
of radius in an accreting source, then we may expect the flow to approximate the 
self-similar solutions discussed here (cf. Fig. 1 in NY). The various properties of 
these solutions then become relevant. To summarize: 

1. The flow is definitely not disk-like in morphology. In fact, the closest analog 
to our solutions in the accretion literature is Bondi (1952) spherical accretion. 
However, our flows differ in important ways from the Bondi problem. The gas 
in our solutions does rotate, it has non-trivial viscous interactions through which 
angular momentum is transported outward, and it has energy dissipation just as 
in ordinary accretion disks. 

2. The angular velocity is significantly sub-Kepler ian and this may have important 
implications for the spin-up of accreting stars. Stars which spin up through an 
advection-dominated mode of accretion are likely to reach a steady state with a 
rotation rate much below the "break-up limit." This may be one solution to the 
angular momentum problem which is discussed frequently in the context of star 
formation. 

3. The radial accretion velocity is typically high in advection-dominated flows. 
Roughly, the velocity scales as f ~ aCs- In a cooling-dominated disk, c^ <^ t;//, 
where Vff is the free-fall velocity, and so v <^ Vff. However, advection-dominated 
flows have Cg ~ Vff and so for a reasonable a ~ 0.1, we find v ~ 0.1^;//. 

4. The Bernoulli parameter is positive over most of the flow (except the regions 
very close to the inner edge), and there is also likely to be violent convection 
especially close to the rotation axis. A fairly substantial bipolar outflow is quite 
likely under these conditions provided only that the material can remain adiabatic 
during the acceleration phase of the outflow. 

5. The convective motions will transport both energy and angular momentum 
and will be a source of viscosity. This has the advantage that, even if for some 
reason the Balbus & Hawley (1991) instability were to be ineffective, convective 
viscosity can still keep the accretion going. 

6. By definition, advection-dominated accretion systems are under-luminous rel- 
ative to the mass accretion rate. This is because the energy is carried along with 
the accreting gas as heat instead of being radiated. If the accreting object is a 
regular star, this energy must finally be radiated from the star and will be seen 
as stellar emission. However, if the accretor is a black hole, then most of the ac- 
cretion energy can disappear through the horizon. It would be very misleading to 
estimate the mass accretion rate of such a system from the observed luminosity. 
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7. The spectrum of an advection-dominated flow is likely to be quite different from 
that of a cooling-dominated accretion disk. If the flow is in the low M optically 
thin limit, then the temperature of the emitted radiation will be close to virial 
and the spectrum will be unusually hard. This appears to be the situation in 
cataclysmic variable boundary layers for M < lO~^'^M0yr~^ (Narayan & Popham 

1993). In the opposite radiation-trapped limit at high M, the photosphere will 
be much farther out than usual and the spectrum will be unusually soft. This is 
apparently the case in FU Orionis systems (Popham et al. 1993). Similar effects 
must be present in accreting neutron star and black hole systems, but this remains 
to be investigated. 

We conclude with a final speculative comment on accretion disk coronae. 
An interesting feature of our solutions in the limit of large e' (efficient cooling) 
is that they have low density hot gas on top of an equatorial thin disk (see 
the e' = 10 solution in Fig. 1). The hot gas resembles the ad hoc coronae that 
various researchers have invoked in models of accretion disks. We speculate (i) 
that coronae are a natural and inevitable feature of any thin disk, (ii) that such 
coronae are best described as advection-dominated fiows rather than as static 
atmospheres, and (iii) that these coronae will themselves have outfiows for the 
reasons discussed in this paper, though with much smaller mass loss rates than in 
fully advection-dominated fiows. Whenever a disk plus corona structure is formed, 
the accreting material must be divided in some self-consistent manner between 
the thin cool disk and the thick hot corona. One possibility is that the corona 
always has just enough mass in it to be marginally advection-dominated, i.e. to 
have f corona ~ 1/2. Detailed calculations with radiative transfer are needed to 
confirm whether such a model will be self-consistent. 
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APPENDIX A: Comparison with Begelman and Meier (1982) 

Begelman & Meier (1982, hereafter BM) considered a problem very similar 
to the one described in this paper, and it is useful to compare the two approaches. 
BM restricted their attention to a radiation-dominated accretion flow which is 
fully advection-dominated. In the language of this paper, their model corresponds 
to the particular case, 7 = 4/3, e' = e = 1. They considered however a more 
general viscosity model than we did, viz. 

uocp'^'p^'r^'gie), (Al) 

where a', /?', 7' are three arbitrary exponents. g{9) is a general function of 
9, though in actual practice, BM set g{9) = constant in their calculations. By 
imposing additional requirements like self-similarity, BM restricted their viscosity 
model to satisfy 

/3' = -a\ 1 + a' - 7' = 1/2, (A2) 

so that they had effectively a one-parameter family of viscosity models parametrized| 
by a'. Our model of viscosity, eq. (2.6), corresponds to the particular choice, 

a' = 1, (3' = -1, 7' = 3/2, {A3) 

and is one member of this family. We have also carried out a few calculations 
with another model corresponding to a' = 1/2, /3' = —1/2, 7' = 1. 

Starting with the basic equations (2.7) - (2.10) which we have written in 
sec. 2, BM assumed that 

Vr, Vq < V^, {AA) 

and derived the following second-order differential equation for O(^), 
(1 - sin^ 9n^) (n" + 3 cot 9Q' - -^ = ^^^^ (%^ + O^^ 

-2{/3' + 1) sin^ ^0'2 _ ^y ^ 3^/ ^ 4^ ^Qg g gj^ 9n'^n'. (Ab) 

The parameter // is the radial exponent of the pressure, p (x r^. For self-similar 
solutions, n = —5/2. Notice that eq. (A5) is only a second-order differential 
equation, whereas the system of equations we derived in this paper is sixth-order. 
The simplification arises because of the assumption Vr ^ t'<^ (eq. A4). 

It is straightforward to derive an equation similar to (A5) starting with our 
equations (2.16) - (2.19). As discussed in the paper, our solutions have v ^^ a. 
Therefore the condition v <^ v^j) is equivalent to the assumption a <^ 1. In fact, 
all we need is a^ <^ 1. Under this condition, we can neglect all terms involving 
v"^ and av in eqs. (2.16), (2.17) and (2.19). Eq. (2.16) then gives 

cl = -{l-siu^9Q^), (A6) 

5 
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while (2.17) gives 



A (pc2) = cos e sin Opn'^. (A7) 

du 



Eliminating v between (2.18) and (2.19) and substituting (A6) and (A7) we then 
find 

(1 - sin^ eif){9," + 3 cot Oil' -\^) = — ^ sin^ Oil (j9? + 9,'A 

-- cos ^ sin ^0^0'. (A8) 

This is exactly equivalent to eq. (A5) provided we set e' = 1 (radiation-dominated, 
full advection) and (3' = —1, 7' = 3/2 (eq A3), n = —5/2. We thus find perfect 
agreement between our equations and those of BM so long as a^, v^ <^1. 

All the solutions we have presented in this paper have a = 0.1 or smaller, 
so that a^ <^ 1. Furthermore, all of our solutions have well-behaved non-singular 
v{9) and satisfy t;^ ^ 1 at all 9. Therefore, although our solutions were obtained 
by solving the exact sixth-order equations (2.16) - (2.19), they do in fact satisfy 
BM's second-order equation (A5) very accurately. These solutions seem to have 
been missed by BM. 

The solutions that BM described in their paper correspond to fiows which 
extend from the equatorial plane at 6* = 7r/2 to a free surface at 9 = 9f, inside 
of which is an empty funnel. A somewhat disturbing feature of their solutions 
is that the radial velocity v diverges at 9 = 9f- The divergence arises because 
the ram pressure term v'^ /2 and the poloidal viscous term {d/d9){adv/d9) in eq. 
(2.16), which would normally control the divergence, were eliminated through the 
assumption Vr ^ V(f,. A divergent solution is of course not compatible with the 
original assumption v <^ 1, and this implies an inconsistency in the solutions. 
Another odd feature of the BM solutions is that for fixed values of the parameters 
(e', a, a', /3', 7', /u), they find a continuous family of self-similar solutions, with 
a continuously tunable funnel opening angle 9f- We find this infinity of solutions 
somewhat disturbing and worry that perhaps a boundary condition at the free 
surface may have been missed. We note that in our calculations, when we fix a 
and e', we find a unique fundamental solution. (We do have the curious higher- 
order solutions described in §3.1, where Q{9) has one or more nodes, but these 
solutions still form only a discrete set, not the continuous family that BM find.) 

We feel that it is important to search for solutions with empty funnels and free 
surfaces using the full sixth-order set of equations (2.16) - (2.19). The numerical 
methods which we have employed in the calculations described in this paper are 
incapable of finding such solutions. It is therefore a completely open question 
whether or not such solutions exist. Tori and empty funnels have been much 
discussed in accretion astrophysics, expecially in the context of AGN (see Frank, 
King & Raine 1992 for a review). However, except for the work of BM, there have 
been no self-consistent dynamical models of such fiows which include viscosity 
fully. In our opinion this is an important problem. 
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APPENDIX B: Non- Rotating Spherical Inflow 

In addition to tiie rotating solutions wtiicti are the main focus of the paper, 
the self-similar equations (2.16)-(2.19) allow a second non-rotating spherically 
symmetric branch of solutions which is closely related to the Bondi (1952) prob- 
lem. To derive this solution we set O = and d/d6 = 0. Equations (2.17) and 
(2.18) then immediately drop out. Equation (2.19) gives 

-— = 3„^ (Bi) 

which has two solutions for v. One of these is t; = 0, which is the condition 
satisfied at 6* = by all the rotating solutions described in the paper. For spherical 
infiow, we consider the second possibility, namely 

Substituting this in eq (2.16), we then solve for the sound speed, to obtain 

Equations (B2), (B3) coupled with the self-similar scalings completely describe 
this branch of solutions. We see that these solutions are allowed whenever |e'| < 

2\/2a (or e' < 5 if a is large). Thus, a self-similar form of spherical accretion is 
possible so long as I7 — 5/3| < a. 

Note that in pure Bondi (1952) spherical fiow, a self-similar form of accretion 
or outflow is allowed only for a single value of 7, viz. 7 = 5/3. In contrast, we 
find that self-similarity is possible in the present problem for a range of values 
of 7. The additional freedom arises because we have viscosity. Even though the 
fiow is not rotating, it still has a velocity divergence and this gives rise to the 
viscous dissipation term ?>v'^ in the right of eq (Bl). When 7 7^ 5/3, self-similarity 
requires a radial gradient in the entropy (as discussed in §3.1). The solution feeds 
this gradient by tuning the magnitude of v so as to supply the required energy 
input through viscous dissipation. 

The other interesting feature is that we can have either self-similar spherical 
infiow or outfiow, depending on the sign of e' (eq B2). Since viscosity is always 
a source of energy, the entropy has to increase in the direction of the fiow. 
For 7 < 5/3, the entropy increases inward in the self-similar solution, and so the 
motion is also inwards, i.e. we have accretion. However, for 7 > 5/3, the solution 
corresponds to entropy increasing outward. In this case, therefore, our solution 
corresponds to a spherical self-similar wind. 

Note that the spherical solutions described here may be either subsonic or 
supersonic. For e' close to zero, t; <^ Cg, and we have subsonic conditions. How- 
ever, as |e'| — s> 2\/2a, the radial velocity increases and the sound speed decreases 
and we can have supersonic fiows. 

As a final comment, we recall that the rotating solutions described in the 
main text of the paper exist only for t' > Co? where C is a constant of order 
unity (§3.1). On the other hand, the spherical solutions described here exist for 

e' < 2^/2a. Thus, for Cc^ < e' < 2\f2a both solution branches are allowed. What 
is the relationship between the two solutions? Which solution will an actual fiow 
prefer? Perhaps the answer is determined by the initial conditions of the fiow, 
especially the angular momentum of the accreting gas. 

22 



References 

Abramowicz, M., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 
646 

Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214 

Begelman, M. C. 1978, MNRAS, 184, 53 

Begelman, M. C. & Meier, D. L. 1982, ApJ, 253, 873 

Bondi, H. 1952, MNRAS, 112, 195 

DubruUe, B. 1992, A&A, 266, 592 

Eggum, G. E., Coroniti, F. V., & Katz, J. I. 1988, ApJ, 330, 142 

Frank, J., King, A., & Raine, D. 1992, Accretion Power in Astrophysics (Cam- 
bridge, UK: Cambridge Unniversity Press) 

Hawley, J. F., Gammie, C. F. & Balbus, S. A. 1994, preprint 

Henriksen, R. N. & Valls-Gaboud, D. 1994, MNRAS, 266, 681 

Landau, L. D. & Lifshitz, E. M. 1959, Fluid Mechanics (London: Pergamon Press) 

Liang, E. P. 1988, ApJ, 334, 339 

Lynden-BeU, D. & Pringle, J. E. 1974, MNRAS, 168, 603 

Mihalas, D. & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics 
(New York: Oxford University Press) 

Narayan, R., Loeb, A., & Kumar, P. 1994, ApJ, 431, 359 

Narayan, R. & Popham, R. 1993, Nature, 362, 820 

Narayan, R. & Yi, I, 1994, ApJ, 428, L13 

Popham, R., Narayan, R., Hartmann, L., & Kenyon, S. 1993, ApJ, 415, L127 

Press, W. H., Teukolsky, S. A., Vetterhng, W. T., & Flannnery, B. P. 1992, 
Numerical Recipes (Cambridge: Cambridge University Press) 

Rees, M. J., Begelman, M. C, Blandford, R. D., & Phinney, E. S. 1982, Nature, 
295, 17 

Ryu, D., & Goodman, J. 1992, ApJ, 388, 438 

Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337 

Spruit, H. C, Matsuda, T., Inoue, M., & Sawada, K. 1987, MNRAS, 229, 517 



23 



Figure Captions 

Figure 1. Self-similar solutions corresponding to a = 0.1, e' = 0.1, 1, 10. Top left: 
angular velocity O as a function of polar angle 9. Top right: radial velocity v. 
Bottom left: density p. Bottom right: square of the sound speed, c^. 

Figure 2. Isodensity contours in the meridional plane for four solutions. The top 
two panels and the bottom left panel correspond to the solutions shown in Figure 
1. The bottom right panel shows a solution in which O reverses sign twice (see 
text). 

Figure 3. Comparison of the exact solutions of this paper with the height- inte- 
grated solutions of NY. The solid lines correspond to the height-integrated solu- 
tions; the dotted lines correspond to spherical averages (cf eq. 3.5); the dashed 
lines correspond to cylindrical ^-averages (eqs. 3.6-3.8); and the long-dashed lines 
correspond to midplane values. The height-integrated values agree remarkably 
well with the spherical averages for all values of e'. 

Figure 4. Dimensionless Bernoulli parameter 6 as a function of the spherical polar 
angle 9 for self-similar solutions with a = 0.1, e = 0.333 (i.e. 7 = 1.5), and from 
bottom to top, e' = 300, 100, 10, 1, 0.75, 0.33 (i.e. / = 0.0011, 0.0033, 0.033, 
0.33, 0.4465, 1). Note that the Bernoulli parameter is always positive close to 
the rotation axis, and is positive at all 9 for e' < 0.75 (i.e. / > 0.446). 

Figure 5. The ratio between the convection term and advection term in the energy 
equation, shown as a function of 6*, for a solution with a = 0.01, ac = 0.005, 
e = 0.01, / = 1. For all 9 > 9crit = 33°, the convection term is smaller than the 
advection term and self-similar advection-dominated accretion is possible. For 
9 < 9crit, convection dominates over advection and we speculate that convection 
will initiate a bipolar outflow in this region of the flow. Note from Fig. 4 that this 
region has the most positive Bernoulli parameter and is therefore most susceptible 
to being ejected. 
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